$ontext
This is the baseline model for dissertation modeling.
$offtext

$set scenario E8_con_50C_SC2

Option   Solprint = off;
Option   Limrow = 0;
Option   Limcol = 0;

sets      plot   biomass source locations  /
$include wsc_pseudocounties_E8_con.csv
/
         brfn    biorefinery sites /
$include brfn_all.csv
/
         feed    /miscanthus, switchgrass/
         crop /ncrp, soy, srg, wht, pnt, ric, osg, cot, hay, oat, crn, brl, bns, sbt, sfl/
         crop_cat /ncrp, crp/
         extra   feedstock classification        /extra1, extra2/
         comp    feedstock composition component /cellulose, hemicellulose, lignin, HHV/
         tech         tech type   /lce, ft_jet/
         term            terminal locations  /
$include terminal_list.csv
/
         trans        transportation mode /hours, dist, total_cost/;

set      state           /
$include state_list.csv
/;


alias (tech, tech2);

Parameters       pbcost(plot, brfn)        '$ per ton for feedstock transport'  /
$ondelim
$include regional_links_E8_con.csv
$offdelim
/

                 btcost(brfn, term)

         acres(plot, crop) 'matrix of supply quantities' /
$ondelim
$include wsc_acres_E8_con.csv
$offdelim
/
         wta(plot, feed, crop)    'values for price levels'    /
$ondelim
$include wsc_wta_E8_con.csv
$offdelim
/
         yield(plot, feed)       /
$ondelim
$include wsc_yields_E8_con.csv
$offdelim
/

         feed_MC(feed)                   moisture content of dry biomass feedstocks   /
miscanthus           0.15
switchgrass              0.15/

         harvest_loss            fraction of total biomass yield not harvested    /0.1/
         logistics_loss          fraction of biomass lost in logistics           /0.05/

         conv_eff(feed, tech)            conversion efficiency for converting feedstock to fuel

         conv_factor(tech, comp)    converts conversion efficiencies into gallons product per ton feedstock             /
lce.hemicellulose    176.9
lce.cellulose        172.85
ft_jet.HHV        7.102/
         tech_eff(tech, comp)       conversion efficiency for cellulosic technologies     /
lce.hemicellulose    0.765
lce.cellulose        0.72
ft_jet.HHV        0.465/

         electricity(tech, feed)               electricity produced per ton of biomass
         elec_eff(tech)          /
lce              0.0366
ft_jet        0.0032/
         naptha(tech, feed)                    naptha produced per ton of biomass
         naptha_eff(tech)                /
ft_jet        0.0/
         gge_conversion(tech)          converts fuel quanitities into equivalent gallons of gasoline   /
lce              0.657
ft_jet        1/

         vmt_fraction(term)      'fraction of 2015 vmt allocated to each terminal'  /
$ondelim
$include terminal_vmt_fraction.csv
$offdelim
/
         fuel_price_base              fuel price for each model run                                  /2.8/
         fuel_price(term)        terminal specific fuel price accounting for historic price differentials
         naptha_price
         elec_price              price of electricity in dollars per kWh                         /0.05/
         carbon_price            price of carbon in dollars per metric ton CO2e                  /50/
         total_supply            target supply quantity for cost min model                      /3000/
         marg_cost_land(crop)
         distance(plot,brfn)
;

naptha_price = 0.86*fuel_price_base;
*adder for additional land rent on cropland which is $5 per thousand acres
marg_cost_land(crop) = 5/1000/10000;
marg_cost_land('ncrp') = 0;
yield(plot, feed) = yield(plot,feed)*(1 - harvest_loss);
wta(plot,feed,crop) = wta(plot,feed,crop)/(1 - harvest_loss);

Table    biomass_composition(feed, comp)
$ondelim
$include feedstock_composition.csv
$offdelim
;

Table    feedstock_od_table(plot, brfn, trans)
$ondelim
$include p2b_limited_con.txt
$offdelim
;

Table        terminal_od_table(brfn, term, trans)
$ondelim
$include wsc_brfn2term.csv
$offdelim
;
Table            terminal_cost(term, tech)
$ondelim
$include terminal_cost.csv
$offdelim
;

Parameter state_terminals(term, state)      /
$ondelim
$include terminal_fips.csv
$offdelim
/
         state_brfn(brfn, state)      /
$ondelim
$include brfn_fips.csv
$offdelim
/
         state_price_delta(state)        /
$ondelim
$include state_price_delta.csv
$offdelim
/
         crp_cat_matching(crop_cat, crop)        /
$ondelim
$include crop_cat_crosswalk.csv
$offdelim
/
         regional_locations(brfn, tech, extra)   /
$ondelim
$include locations_E8_con_r1a.put
$include locations_E8_con_r1b.put
$include locations_E8_con_r2.put
$include locations_E8_con_r3.put
$include locations_E8_con_r4.put
$include locations_E8_con_r5.put
$include locations_E8_con_r6a.put
$include locations_E8_con_r6b.put
$include locations_E8_con_r7.put
$include locations_E8_con_r8.put
$include locations_E8_con_r9a.put
$include locations_E8_con_r9b.put
$offdelim
/
;

*pbcost(plot, brfn) = feedstock_od_table(plot, brfn, 'total_cost');
btcost(brfn, term) =  terminal_od_table(brfn, term, 'total_cost');
*btcost(brfn, term)$(btcost(brfn, term) lt 0.5) = 500;
distance(plot,brfn) = feedstock_OD_table(plot,brfn,'dist');

fuel_price(term) = fuel_price_base;
* + sum(state, state_terminals(term, state)*state_price_delta(state));



conv_eff(feed, tech) = sum(comp, biomass_composition(feed, comp)*tech_eff(tech, comp)*conv_factor(tech, comp))/1000;
electricity(tech, feed) = elec_eff(tech)*biomass_composition(feed, 'HHV')/1.1023/.0036/1000;
naptha(tech, feed) = naptha_eff(tech)*biomass_composition(feed, 'HHV')/1.1023/0.1216089/1000;

*Carbon accounting parameters
Parameters
* LUC carbon from GREET
         LUC_CO2(plot, feed, crop)
         LUC_by_cat(plot, feed, crop_cat)    'LUC emissions from GREET Mg C/acre'    /
$ondelim
$include C_dLUC.csv
$offdelim
/
         cultivation_co2(plot, feed)       'biomass cultivation emissions kg CO2e/dry ton sold'
         c02_per_ton_mile        'Carbon emissions in biomass transport Mg CO2e/dry ton-mile'  /0.000177/
         conversion_co2(tech)    'Carbon emissions at brfn Mg CO2e/dry ton consumed'    /
lce              0
ft_jet        0.01775/
         CO2_elect(brfn, tech)   'Carbon emissions of displaced electricity kg/MWh'
         state_elect_CI(state)   'Carbon intensity of electricity in each state in kg/MWh from 2019 eGRID' /
$ondelim
$include state_electricity_CI.csv
$offdelim
/
         fuel_transport_CO2      'Carbon emissions kg per 1,000 gallon-miles of fuel transport' /0.645/
         bt_CO2(brfn,term)       'fuel transport emissions kg CO2e/1,000 gal-miles'
         end_use_CO2(tech)       'Carbon emissions at the end use kg/1,000 gallons'  /
lce              0.0366
ft_jet        0.0032/
;

* Calculate emissions over each fuel link
bt_CO2(brfn, term) =  fuel_transport_CO2*terminal_od_table(brfn, term, 'dist');

* Assign state electricity CI to brfns in the state, convert to Mg/Mwh
CO2_elect(brfn, tech) = (sum(state, state_brfn(brfn, state)*state_elect_CI(state)))/1000;

* Calculate cultivation emissions based on yield from crop budgets and GREET1_2021, convert to Mg/ton
cultivation_co2(plot, 'switchgrass')$(yield(plot, 'switchgrass') gt 0.01) = (38.56 + 13.56/yield(plot, 'switchgrass'))/1000;
cultivation_co2(plot, 'miscanthus')$(yield(plot, 'miscanthus') gt 0.01) = (28.37 + 2.97/yield(plot, 'miscanthus'))/1000;

*Convert soil organic carbon changes into CO2 emissions
LUC_CO2(plot, feed, crop) = - sum(crop_cat, crp_cat_matching(crop_cat, crop)*LUC_by_cat(plot, feed, crop_cat))*44/12;

*** CHECK UNITS TO GET CONSISTENT and APPROPRIATE CARBON VALUES

*Removing unneeded varibales and equations CHECK BEFORE RUNNING!
sets spar(plot, feed, brfn)  set used to eliminate unnecessary links
     spar1(plot, feed, crop)       set used to eliminate unnecessary supplies
         spar2(brfn, tech, extra)
         spar4(feed, brfn, tech, extra)
         spar3(brfn, tech)
         spar5(brfn, term, tech)
         spar6(feed, brfn)
         spar7(plot, crop)
         spar8(plot, feed);

spar(plot, feed, brfn)$(pbcost(plot, brfn)*sum(crop, acres(plot, crop))*yield(plot, feed) > 0.5)=yes;
spar1(plot, feed,  crop)$(acres(plot, crop)*sum(brfn, pbcost(plot, brfn))*yield(plot, feed) > 0.5)=yes;
spar2(brfn, tech, extra)$(regional_locations(brfn, tech, extra) gt 0.5) = yes;
spar3(brfn, tech)$(sum(extra, regional_locations(brfn, tech, extra)) gt 0.5) = yes;
spar4(feed, brfn, tech, extra)$(regional_locations(brfn, tech, extra) gt 0.5) = yes;
spar5(brfn, term, tech)$(sum(extra, regional_locations(brfn, tech, extra)) gt 0.5) = yes;
spar6(feed, brfn)$(sum((tech, extra), regional_locations(brfn, tech, extra)) gt 0.5) = yes;
spar7(plot, crop)$(acres(plot, crop)*sum(brfn, pbcost(plot, brfn))*sum(feed, yield(plot, feed)) gt 0.5)=yes;
spar8(plot, feed)$(sum(crop, acres(plot, crop))*sum(brfn, pbcost(plot, brfn))*yield(plot, feed) > 0.5)=yes;


*Assign missing links a very high cost
pbcost(plot, brfn)$(pbcost(plot, brfn) lt 0.05) = 200;


*Fix the parameters to have units of $10M, 10M gallons, 10k tons and 10k Mg Co2e
pbcost(plot, brfn) = pbcost(plot, brfn)/1000;
acres(plot, crop) = acres(plot, crop)/10000;
btcost(brfn, term) = btcost(brfn, term)/100;
wta(plot, feed, crop) = wta(plot, feed, crop)/1000;
carbon_price = carbon_price/1000;

Positive Variables
         adopted_acres(plot, feed, crop)         area planted in energy crop "feed" in "plot" with incumbent "crop"
         biomass_consumed(feed, plot, crop)      annual quantity of feedstock consumed from plot "i" at price level "plev"
         p2b(feed, plot, brfn)                   annual quantity of biomass from plot "i" taken to biorefinery "k"
         Qb(brfn, tech, extra)                   quantity of fuel "b" produced at biofinery "k"
         T(brfn, term, tech)                     quantity of fuel delivered to terminal from biorefinery
         x1(feed, brfn, tech, extra)             linearizing variables for biorefinery costs
         crop_displaced(crop)                    total displaced area for a given crop
*         total_supply
;

Variable
         total_cost
         total_carbon
         profit;

binary variable
         xi(brfn, tech, extra);

Equations
         annual_cost                     total annual cost for biofuel supply
         annual_profit
         xi_constraint(brfn, tech, extra)       constraint for integer variable xi
         Qb_constraint(brfn, tech, extra)       relates Qb to x1's

*Constraints to set up biofuel networks
         land_constraint(plot, crop)             limits conversion to available land and prevents double counting of land
         supply_constraint(plot, feed, crop)     limits biomass leaving a plot at a price level to what is available
         plot_flow(feed, plot)                   constrains feedstock leaving plot to what is consumed
         brfn_cap(feed, brfn)                   limits fuel production capcity at biorefinery to be less than or equal to the feedstock inputs
         brfn_flow(brfn, tech)                   limits the fuel leaving the brfn to be less than the fuel produced
*         total_demand                            aggregates the fuel quantities to the set limit
*         terminal_min(term, fuel)                lower bound on fuel delivered to terminal
*         terminal_max(term)                upper bound on fuel delivered to terminal
*         terminal_max2(term)
*         terminal_max3(term)
*         feed_utilization(feed)                  upper bound on feedstock use nationally
         crop_accounting(crop)
         carbon_accounting

*Redundant constraints that may better define the model
*         redundant1(brfn, extra)
         redundant2(brfn, tech)
;


Parameters
         ac(tech)                         fixed cost component in linearized biorefinery costs      /
lce              3.6917481
ft_jet        33.1629/
         bc(tech)                         feed capacity depedent component in linearized biorefinery costs      /
lce              0.4528
ft_jet        0.48858/
         iir             /0.1/
         lifetime        /20/
         CRF
         ao(tech)                         fixed cost component in linearized biorefinery costs      /
lce              0.4356352
ft_jet        3.15/
         bo(tech)                         tech cap dependent component in linearized biorefinery costs    /
lce              0.05952
ft_jet        0.0764/
         MM(tech)                               Maximum capacity                                  /
lce              131
ft_jet        200/
         a(brfn, tech, extra)
         b(brfn, tech, extra)
         M(brfn, tech, extra);

CRF = iir*(1+iir)**lifetime/((1+iir)**lifetime - 1);


a(brfn, tech, extra) = (CRF*ac(tech) + ao(tech));
a(brfn, tech, 'extra2') = (CRF*ac(tech) + ao(tech)) + 0.1;
b(brfn, tech, extra) = (CRF*bc(tech) + bo(tech));
M(brfn, tech, extra) = MM(tech);

*Objective
annual_cost..            total_cost =e= sum((spar(plot, feed, brfn)), pbcost(plot, brfn)*p2b(feed, plot, brfn)) + sum((spar1(plot, feed, crop)), wta(plot, feed, crop)*biomass_consumed(feed, plot, crop))
                         + sum((feed, brfn, tech, extra), b(brfn, tech, extra)*x1(feed, brfn, tech, extra)) + sum((brfn, tech, extra), a(brfn, tech, extra)*xi(brfn, tech, extra))
                         + sum((brfn, term, tech), T(brfn, term, tech)*(btcost(brfn,term)+ terminal_cost(term, tech))) + sum(crop, marg_cost_land(crop)*crop_displaced(crop));

annual_profit..          profit*1000 =e= sum((brfn, term, tech), fuel_price(term)*gge_conversion(tech)*T(brfn, term, tech))+ elec_price*sum((brfn, tech, extra, feed), electricity(tech, feed)*x1(feed, brfn, tech, extra))
                         + naptha_price*sum((brfn, tech, extra, feed), naptha(tech, feed)*x1(feed, brfn, tech, extra)) - total_cost - carbon_price*total_carbon/1000;

*Subject to:
land_constraint(plot, crop)$spar7(plot, crop)..    sum(feed, adopted_acres(plot, feed, crop)) =l= acres(plot, crop);
supply_constraint(plot, feed, crop)$spar1(plot, feed, crop)..     biomass_consumed(feed, plot, crop) =e= yield(plot, feed)*adopted_acres(plot, feed, crop);
plot_flow(feed, plot)..        sum(spar1(plot, feed, crop), biomass_consumed(feed, plot, crop)) =g= sum(spar(plot, feed, brfn), p2b(feed, plot, brfn));
brfn_cap(feed, brfn)..        sum(spar(plot, feed, brfn), p2b(feed, plot, brfn)*(1-logistics_loss))
                                                                 =g= sum((tech, extra), x1(feed, brfn, tech, extra));
brfn_flow(brfn, tech)..          sum(term, T(brfn, term, tech)) =l= sum(extra, Qb(brfn, tech, extra));
*terminal_max(term)..       sum(brfn, T(brfn, term, 'lce')) =l= vmt_fraction(term)*Ethanol_demand;
*terminal_max2(term)..       sum(brfn, T(brfn, term, 'ft_diesel')) =l= vmt_fraction(term)*Diesel_demand*0.95;
*terminal_max3(term)..       sum((brfn, tech), T(brfn, term, tech)) =l= vmt_fraction(term)*3300 ;
crop_accounting(crop)..            sum(spar8(plot, feed), adopted_acres(plot, feed, crop)) =e= crop_displaced(crop);
*feed_utilization(feed)..         sum(spar1(r_plot, feed, plev), biomass_consumed(feed, r_plot, plev)) =l= resource_limit(feed)*sum((r_plot,plev), supply(r_plot, feed, plev));

*Demand constraint
*total_demand..           sum((brfn, term, tech), gge_conversion(tech)*T(brfn, term, tech)) =g= total_supply;

xi_constraint(brfn, tech, extra)..      sum(feed, x1(feed, brfn, tech, extra)) =l= M(brfn, tech, extra)*xi(brfn, tech, extra);
*redundant1(brfn, extra)..         sum(tech, xi(brfn, tech, extra)) =l= 1;
redundant2(brfn, tech)..         sum(term, T(brfn, term, tech)) =l= sum(extra, M(brfn, tech, extra)*xi(brfn, tech, extra));
Qb_constraint(brfn, tech, extra)..      Qb(brfn, tech, extra) =l= sum(feed, x1(feed, brfn, tech, extra)*conv_eff(feed, tech));

carbon_accounting..      total_carbon =e= sum((spar1(plot, feed, crop)), LUC_CO2(plot, feed, crop)*adopted_acres(plot, feed, crop) + cultivation_co2(plot, feed)* biomass_consumed(feed, plot, crop))
                         + sum(spar(plot, feed, brfn), p2b(feed, plot, brfn)*distance(plot, brfn)*c02_per_ton_mile)
                         + sum(spar4(feed, brfn, tech, extra), conversion_co2(tech)*x1(feed, brfn, tech, extra)) - sum(spar4(feed, brfn, tech, extra), CO2_elect(brfn, tech)*electricity(tech, feed)*x1(feed, brfn, tech, extra))
                         + sum(spar5(brfn, term, tech), T(brfn, term, tech)*(bt_CO2(brfn,term)+ end_use_CO2(tech)));

xi.up(brfn, tech, extra) = regional_locations(brfn, tech, extra);
*T.up(brfn, term, 'ft_diesel') = 0;
*x1.up('msw_dirty', brfn, 'ft_diesel') = 500;



*
*xi.up(brfn, tech, extra) = regional_locations(brfn, tech, extra);
adopted_acres.up(plot, feed, crop) = acres(plot, crop);
*T.up(brfn, term, 'ft_jet') = 0;
*x1.up('msw_dirty', brfn, 'ft_jet') = 500;

*xi.up(brfn, 'lce', extra) = 0;

Option MIP = cplex;
Option iterlim = 1000000;
Option reslim = 18000;

*Execute_loadpoint 'USDA_p';

Model USDA /ALL/;
*$ontext
FILE opt "Cplex option file" / cplex.opt /;
PUT opt;
PUT
'workmem 24000'/
'memoryemphasis 1'/
*'nodefileind '/
*'nodesel     2'/
*'varsel      3'/
'brdir 1'/
'cuts    2'/
*'cutlo 0.0123'/
*'mipemphasis 3'/
*'mipstart 1'/
'probe   3'/
*'probetime 3600'/
*'indic redundant1(brfn, tech)$xi(brfn, tech) 0'/
'parallelmode -1'/
'threads=0'/;
PUTCLOSE OPT;

USDA.optfile=1;
USDA.OptCR = 0.015;
*$offtext
*Option savepoint=1;

*$ontext
Set      run  /run0*run20/ ;

alias  (extra,extras);
alias    (term, term2);

Parameter fprice(run)   defines the fraction of the max demand for each run /
run0        3.8
run1        3.9
run2        4
run3        4.1
run4        4.2
run5        4.3
run6        4.4
run7        4.5
*$ontext
run8        4.6
run9        4.7
run10        4.8
run11        4.9
run12        5
run13        5.1
run14        5.2
run15        5.3
run16        5.4
run17        5.5
run18        5.75
run19        6
run20        7
*$offtext
/;


*$ontext
*Create results output files
file status_%scenario%;
file results_%scenario%   ;
file results_%scenario%_brfn;
*file results_%scenario%_feedstock_links;
*file results_%scenario%_fuel_links;
file results_%scenario%_pc;
*file update;

*Make results files CSV
results_%scenario%.pc=5;
results_%scenario%_brfn.pc=5;
*results_%scenario%_feedstock_links.pc=5;
*results_%scenario%_pc.pc=5;
results_%scenario%_brfn.pw= 1500;
results_%scenario%.pw=1500;

$onend


$ontext
put results_%scenario%_feedstock_links;
put "scenario", "fuel_price ($/gge)", "source_id", "dest_id", "energy_crop", "quant_tons"/;

put results_%scenario%_fuel_links;
put "scenario", "fuel_price ($/gge)", "source_id", "dest_id", "pathway", "fuel_deliveries (MGY)"/;
$offtext

put results_%scenario%_pc;
put "scenario", "fuel_price ($/gal)", "source_id", "energy_crop", "avg_procurement", "quantity (bdt/yr)", "area (acres)", "LUC (Mg CO2e/yr)", "cultiv_GHG (Mg CO2e/yr)"/;
put results_%scenario%_brfn;
put "scenario", "fuel_price ($/gal)", "brfn_id", "technology", "extra", "fuel_output (MGY)", "electricity output (GWh/yr)", "naptha (MGY)", "feedstock_cap (bdt/day)",  ;
         loop (feed)     do
         put feed.tl,
         endloop;
put "capital_cost (M$)", "annual_cost (M$/yr)", "annual_capital (M$/yr)", "O&M (M$/yr)", "feedstock_procurement (M$/yr)", "feedstock_transport (M$/yr)", "fuel_distribution (M$/yr)",
         "carbon cost (M$/yr)", "ac ($/gal)", "max_trans_dist (miles)", "totalGHG (Mg/yr)", "CI (g/MJ)"/;

put results_%scenario%;
put "scenario", "fuel_price ($/gal)", "total_biofuel (MGal)", "total_land (Macres)", "total_GHG (mmMg/yr)";
         loop    (crop) do
         put     crop.tl,
         endloop;

         loop    (crop) do
         put     crop.tl,
         endloop;

         loop (tech) do
         put tech.tl,
         endloop;
put /;
$ontext
put " ", "$/gge", "MGGEY", "GWh/yr", "MGY",
         loop (tech, feed)$(feed_tech_match(feed, tech) gt 0.5)        do
         put feed.tl,
         endloop;
         loop (tech)     do
         put "# brfns",
         endloop;
         loop (feed)     do
         put "consumption ktons/yr",
         endloop;
put "billion $", "billion $",;
         loop (tech) do
         put "brfn capital (M$)",
         endloop;
put "million ton miles/yr", "million ton miles/yr", "million ton miles/yr", "million ton miles/yr", "million ton miles/yr)", "million ton miles/yr", "million ton miles/yr", "million ton miles/yr", "million ton miles/yr",
         "million acres", "milion acres"/;
$offtext

Parameter
                 feedstock_consumption(feed)
                 brfn_production(brfn, tech, extra)
                 brfn_consumption(brfn, tech, extra, feed)
*                 links(plot, brfn, feed)
*                 term_links(brfn, term, tech)
                 consumed_biomass(plot, feed, crop)
                 quant
                 biomass_production(crop)
*                 quant_last_run
                 brfn_quant(brfn, tech, extra)
                 brfn_count(tech)
                 elect_quant
                 naptha_quant
                 fuel_by_tech(tech)
                 total_annual_cost
                 brfn_capital_total(tech)
                 brfn_electricity(brfn, tech, extra)
                 brfn_naptha(brfn, tech, extra)
                 feedstock_cap(brfn, tech, extra)
                 capital_cost(brfn, tech, extra)
                 brfn_annual_cost(brfn, tech, extra)
                 annual_capital(brfn, tech, extra)
                 O_M(brfn, tech, extra)
                 feedstock_procurement(brfn, tech, extra)
                 avg_feed_pc(plot, feed)
                 feedstock_transport(brfn, tech, extra)
                 fuel_distribution(brfn, tech, extra)
*                 marginal_feedstock(brfn, tech, extra)
*                 mc(brfn, tech, extra)
                 avg_cost(brfn, tech, extra)
                 max_trans_dist(brfn, tech, extra)
                 crops_displaced(crop)
                 crop_acres_changed(plot, feed, crop)
                 total_land
                 total_GHG
                 brfn_GHG(brfn, tech, extra)
                 brfn_CI(brfn, tech, extra)
                 brfn_carbon_cost(brfn, tech, extra)
                 LUC(plot, feed, crop)
                 cultivation_GHG(plot, feed, crop)
                 avg_feed_GHG(plot, feed)
                 site_biomass(brfn)
                 fff(plot, feed);

$onend

*Solve USDA using MIP maximizing profit;

*quant = 0;

*$ontext
Loop run do

*set new demand constraint
fuel_price_base = fprice(run);
fuel_price(term) = fuel_price_base;
*+ sum(state, state_terminals(term, state)*state_price_delta(state));
naptha_price = fuel_price_base*.86;
*pbcost(plot, brfn, feed) = pbcost_base(plot, brfn, feed) + (fuel_price - 3.55*gge_conversion('fahc'))*diesel_use(plot, brfn, feed)/1000;

*xi.lo(brfn, tech, extra)$(xi.l(brfn, tech, extra) gt 0.75) = 1;

*quant_last_run = quant;

*put update;
*put "Working on " run.tl /
*putclose;

*Resolve for new demand constraint
Solve USDA using MIP maximizing profit;


*$ontext
*Calculate results for export
brfn_production(brfn, tech, extra) = sum(feed, gge_conversion(tech)*conv_eff(feed, tech)*x1.l(feed, brfn, tech, extra))*10;
consumed_biomass(plot, feed, crop) = biomass_consumed.l(feed, plot, crop)*10000;
fff(plot, feed) = sum(crop, consumed_biomass(plot, feed, crop));
site_biomass(brfn) = sum((feed, tech, extra), x1.l(feed, brfn, tech, extra));

loop (plot, feed)$(fff(plot, feed) gt 10) do
         avg_feed_pc(plot, feed) = sum((crop), wta(plot, feed, crop)*biomass_consumed.l(feed, plot, crop) + marg_cost_land(crop)*adopted_acres.l(plot, feed, crop))/sum((crop), biomass_consumed.l(feed, plot, crop));
         LUC(plot, feed, crop) = LUC_CO2(plot, feed, crop)*adopted_acres.l(plot, feed, crop)*10000;
         cultivation_GHG(plot, feed, crop) = cultivation_co2(plot, feed)*consumed_biomass(plot, feed, crop);
         avg_feed_GHG(plot, feed) = (sum(crop, LUC(plot,feed, crop) + cultivation_GHG(plot, feed, crop)))/sum((crop), consumed_biomass(plot, feed, crop));
endloop;

loop (brfn, tech, extra)$(brfn_production(brfn, tech, extra) gt 1) do
         brfn_quant(brfn, tech, extra) = sum(feed, conv_eff(feed, tech)*x1.l(feed, brfn, tech, extra))*10;
         brfn_consumption(brfn, tech, extra, feed) = x1.l(feed, brfn, tech, extra)*10000;
         brfn_electricity(brfn, tech, extra) = sum(feed, electricity(tech, feed)*x1.l(feed, brfn, tech, extra))*10;
         brfn_naptha(brfn, tech, extra) = sum(feed, naptha(tech, feed)*x1.l(feed, brfn, tech, extra))*10;
         feedstock_cap(brfn, tech, extra) = sum(feed, brfn_consumption(brfn, tech, extra, feed))/365/.9;
         capital_cost(brfn, tech, extra) =  (sum(feed, bc(tech)*x1.l(feed, brfn, tech, extra)) + ac(tech)*xi.l(brfn, tech, extra))*10;
         annual_capital(brfn, tech, extra) = CRF*capital_cost(brfn, tech, extra);
         O_M(brfn, tech, extra) = (sum(feed, bo(tech)*x1.l(feed, brfn, tech, extra)) + ao(tech)*xi.l(brfn, tech, extra))*10;

         feedstock_procurement(brfn, tech, extra) = sum((feed)$(x1.l(feed, brfn, tech, extra) gt 0.1), sum((plot)$(p2b.l(feed, plot, brfn) gt 0.001), avg_feed_pc(plot, feed)*p2b.l(feed, plot, brfn)))
                                                         *10*sum(feed, x1.l(feed, brfn, tech, extra)) / site_biomass(brfn);
         feedstock_transport(brfn, tech, extra) = sum((feed)$(x1.l(feed, brfn, tech, extra) gt 0.1), sum((plot), pbcost(plot, brfn)*p2b.l(feed, plot, brfn)))
                                                 *10*sum(feed,x1.l(feed, brfn, tech, extra))/site_biomass(brfn);
         fuel_distribution(brfn, tech, extra) = sum((term), T.l(brfn, term, tech)*(btcost(brfn,term)+ terminal_cost(term, tech)))
                                                 *10*sum(feed,x1.l(feed, brfn, tech, extra))/site_biomass(brfn);
         brfn_annual_cost(brfn, tech, extra) = annual_capital(brfn, tech, extra) + O_M(brfn, tech, extra) + feedstock_procurement(brfn, tech, extra) + feedstock_transport(brfn, tech, extra) + fuel_distribution(brfn, tech, extra);
*         marginal_feedstock(brfn, tech, extra) = smax((plot, feed, crop)$(x1.l(feed, brfn, tech,extra) gt 1 and p2b.l(feed, plot, brfn) gt 0.001 and biomass_consumed.l(feed, plot, crop) gt 0.001), (wta(plot, feed, crop) + pbcost(plot, brfn))/conv_eff(feed, tech)/gge_conversion(tech));
*         mc(brfn, tech, extra) = (annual_capital(brfn, tech, extra) + O_M(brfn, tech, extra))/brfn_production(brfn, tech, extra) + marginal_feedstock(brfn, tech, extra) + smax((term)$(T.l(brfn, term, tech) gt 0.1), (btcost(brfn,term)+ terminal_cost(term, tech)))/gge_conversion(tech);
         max_trans_dist(brfn, tech, extra) = 0;
*smax((plot, feed)$(x1.l(feed, brfn, tech,extra) gt 0.1 and p2b.l(feed, plot, brfn) gt 0.1), feedstock_od_table(plot, brfn, 'dist') );

         brfn_GHG(brfn, tech, extra) = (conversion_co2(tech)*sum(feed,x1.l(feed, brfn, tech, extra)) - CO2_elect(brfn, tech)*sum(feed,electricity(tech, feed)*x1.l(feed, brfn, tech, extra)))*10000;
         brfn_CI(brfn, tech, extra) = (brfn_GHG(brfn,tech,extra) + (sum((feed)$(x1.l(feed, brfn, tech, extra) gt 0.1), sum((plot)$(p2b.l(feed, plot, brfn) gt 0.001), avg_feed_GHG(plot, feed)*p2b.l(feed, plot, brfn)*10000)) +
                                      sum((feed)$(x1.l(feed, brfn, tech, extra) gt 0.1), sum((plot), distance(plot, brfn)*c02_per_ton_mile*p2b.l(feed, plot, brfn)))*10000)*sum(feed,x1.l(feed, brfn, tech, extra))/site_biomass(brfn))
                                         /(brfn_production(brfn, tech, extra)*129.84) ;
         brfn_carbon_cost(brfn, tech, extra) = brfn_production(brfn, tech, extra)*129.84*carbon_price/1000;
         avg_cost(brfn, tech, extra) = (brfn_annual_cost(brfn, tech, extra)- brfn_electricity(brfn, tech, extra)*elec_price + brfn_carbon_cost(brfn, tech, extra))/(brfn_production(brfn, tech, extra) + 0.86*brfn_naptha(brfn, tech, extra));

endloop;

*fuel_production(tech) = sum((feed, brfn, extra)$(x1.l(feed, brfn, tech, extra) gt 0.05), 10*gge_conversion(tech)*conv_eff(feed, tech, extra)*biogenic_energy(feed)*x1.l(feed, brfn, tech, extra));
feedstock_consumption(feed) = sum((brfn, tech, extra), x1.l(feed, brfn, tech, extra))*10000;







*feedstock_consumption(feed) = sum((brfn, tech, extra), x1.l(feed, brfn, tech, extra));
brfn_count(tech) = sum((brfn, extra), xi.l(brfn, tech, extra));
elect_quant = sum((brfn, tech, extra, feed), electricity(tech, feed)*x1.l(feed, brfn, tech, extra))*10;
naptha_quant = sum((brfn, tech, extra, feed), naptha(tech, feed)*x1.l(feed, brfn, tech, extra))*10;
fuel_by_tech(tech) = sum((brfn, extra), brfn_production(brfn, tech, extra) + 0.86*brfn_naptha(brfn, tech, extra));
total_annual_cost = annual_cost.l*10;
brfn_capital_total(tech) =  sum((feed, brfn, extra), bc(tech)*x1.l(feed, brfn, tech, extra)) + sum((brfn, extra), ac(tech)*xi.l(brfn, tech, extra));

*links(plot, brfn, feed) =  p2b.l(feed, plot, brfn)*10000;
*term_links(brfn, term, tech) = T.l(brfn, term, tech)*10;

crops_displaced(crop) = crop_displaced.l(crop)*10000/1000000;
total_land = sum(crop, crops_displaced(crop));
crop_acres_changed(plot, feed, crop) = adopted_acres.l(plot, feed, crop)*10000;
biomass_production(crop) = sum((plot, feed), consumed_biomass(plot, feed, crop)/1000);
total_GHG = total_carbon.l/100;
quant = sum((brfn, tech, extra), brfn_production(brfn, tech, extra) + 0.86*brfn_naptha(brfn, tech, extra));


*Save results
*Model and solver status output file
put status_%scenario%;
put run.tl, "Model status", USDA.modelstat, "Solver status", USDA.solvestat/;

*Summary of national results
put results_%scenario%;
put      "%scenario%", fuel_price_base, quant, total_land, total_GHG,
         loop (crop) do
         put crops_displaced(crop),
         endloop;
         loop (crop) do
         put biomass_production(crop),
         endloop;
         loop (tech) do
         put fuel_by_tech(tech),
         endloop;
put /;

*Site specific supply curve
put results_%scenario%_brfn;
loop     (brfn, tech, extra)$(brfn_production(brfn, tech, extra) gt 1)         do
put "%scenario%", fuel_price_base, brfn.tl, tech.tl, extra.tl, brfn_quant(brfn, tech, extra), brfn_electricity(brfn, tech, extra), brfn_naptha(brfn, tech, extra), feedstock_cap(brfn, tech, extra),
         loop (feed)     do
         put brfn_consumption(brfn, tech, extra, feed),
         endloop
put capital_cost(brfn, tech, extra), brfn_annual_cost(brfn, tech, extra), annual_capital(brfn, tech, extra), O_M(brfn, tech, extra), feedstock_procurement(brfn, tech, extra), feedstock_transport(brfn, tech, extra),
         fuel_distribution(brfn, tech, extra), brfn_carbon_cost(brfn, tech, extra), avg_cost(brfn, tech, extra), max_trans_dist(brfn, tech, extra), brfn_GHG(brfn, tech, extra), brfn_CI(brfn,tech,extra) /
endloop;

$ontext
*Deliveries between plots and biorefinieries
put results_%scenario%_feedstock_links;
loop (plot, brfn, feed)$(links(plot, brfn, feed) gt 100) do
         put "%scenario%", fuel_price_base, plot.tl, brfn.tl, feed.tl, links(plot, brfn, feed)/
         endloop;

put results_%scenario%_fuel_links;
loop (brfn, term, tech)$(term_links(brfn, term, tech) gt 0.1) do
         put "%scenario%", fuel_price_base, brfn.tl, term.tl, tech.tl, term_links(brfn, term, tech) /
         endloop;
$offtext

put results_%scenario%_pc;
loop (plot, feed, crop)$(biomass_consumed.l(feed, plot, crop) gt 0.000100)  do
         put "%scenario%", fuel_price_base, plot.tl, feed.tl, crop.tl, consumed_biomass(plot, feed, crop), crop_acres_changed(plot, feed, crop), LUC(plot, feed, crop), cultivation_GHG(plot, feed,crop)/
         endloop;

if USDA.solvestat = 1 then
         xi.lo(brfn, tech, extra)$(xi.l(brfn, tech, extra) gt 0.5) = 1;
endif;
*xi.lo(brfn, tech, extra)$(xi.l(brfn, tech, extra) gt 0.5) = 1;
*xi.up(brfn, tech, 'extra2')$(xi.l(brfn, tech, 'extra1') gt 0.5) =1;

*xi.lo(brfn, tech, extra)$(xi.l(brfn, tech, extra) gt 0.1) = xi.l(brfn, tech, extra);
*$offtext
endloop;
